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ABSTRACT 

We present a new time-stepping criterion for iV-body simulations that is based on the 
true dynamical time of a particle. This allows us to follow the orbits of particles cor- 
rectly in all environments since it has better adaptivity than previous time-stepping 
criteria used in A'^-body simulations. Furthermore, it requires far fewer force evalua- 
tions in low density regions of the simulation and has no dependence on artificial pa- 
rameters such as, for example, the softening length. This can be orders of magnitude 
faster than conventional ad-hoc methods that employ combinations of acceleration 
and softening and is ideally suited for hard problems, such as obtaining the correct 
dynamics in the very central regions of dark matter haloes. We also derive an eccentric- 
ity correction for a general leapfrog integration scheme that can follow gravitational 
scattering events for orbits with eccentricity e —>■ 1 with high precision. These new ap- 
proaches allow us to study a range of problems in collisionless and collisional dynamics 
from few body problems to cosmological structure formation. We present tests of the 
time-stepping scheme in A'^-body simulations of 2-body orbits with eccentricity e — > 1 
(elliptic and hyperbolic) , equilibrium haloes and a hierarchical cosmological structure 
formation run. 

Key words: methods: iV-body simulations - methods: numerical 



1 INTRODUCTION 



Achieving spatial adaptivity in the evaluation of forces in A'^- 
body simulations is a well-studied problem with many effec- 
tive approaches based on the use of tree structures and mul- 
tipole expansions or nested grids and FFT techniques. Such 
adaptivity in space also comes with a desire to achieve adap- 
tivity in the time integration of these simulations since a 
large dynamic range in density implies a large dynamic range 
in time-scales for self gravitating systems (T ~ 1/VGp ). 
Two very different problems present themselves when try- 
ing to achieve this. First, there are no practical (explicit), 
general purpose (applicable to a wide range of astrophys- 
ical problems), adaptive integration techniques known for 
the A''-body problem which are symplectic. By this we mean 
that the numerical integration is an exact Hamiltonian phase 
flow very close to the phase flow of the continuous system 
under study. This is a very desirable property for following 
systems for longer than a single dynamical time. Such exact 
preservation of the geometrical properties of the dynamical 
system is possible for fixed time-step schemes. For general 
A'^-body simulations with adaptive time-stepping we have 
to resort to approximate symplectic behaviour or preserva- 
tion of time symmetry (a property which is known to lead 
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to very good integration methods). The second problem is 
that of continuously determining the appropriate time-step 
for each particle in the simulation so that the error in the 
integration remains within tolerance while performing the 
fewest possible force evaluations and minimising the compu- 
tational overhead. Resolving this second problem is the main 
focus of this paper and can be considered independently 
from methods symmetrizing the t ime-stepping scheme such 
as presented in lStadell l|200lh and lMakino et al.l l|2006l '). 

There are several known time-step criteria based on dif- 
ferent properties of the simulation (e.g. local density p(r), 
potential ^(r), softening e, acceleration a, jerk a or even the 
velocity v of the particle) tha t are used in state-of-the-art 
numerical codes jQuinn et al. 1997: Stadcl 2001; Aarset^ 
l2003l : iPower et al.H2003l : ISpringelboOSl ). Some of them have 
a physical motivation, others are just a clever combination 
of physical properties in order to obtain a criterion which 
has the physical unit of time. All are an attempt to find an 
inexpensive way of determining an appropriate time-step for 
each particle in the simulation. For cosmological simulations, 
the ad-hoc time-step criterion based on the acceleration and 
the gravitational softening of the particle (T ~ -^e/a ) has 
proven very successful, despite the fact that it is not di- 
rectly related to the dynamical time in these simulations. 
One reason why this time-step criterion is thought to work 
well is that it results in a very tight time-step distribution 
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with very infrequent changes in the time-step of a particle 
in block time-stepping schemes (power of 2 step sizes; see 
also section [22}. This hides the evils due to the first prob- 
lem since the behaviour is more like that of a fixed time- 
stepping scheme than for other criteria. The price, however, 
is that more time-steps are taken in lower density regions 
of the simulation than would seem to be necessary. Further- 
more, while still adequate for the type of simulations that 
have been performed up to now, which adapt the softening 
and mass of particles in the highest density regions and thus 
reduce the time-step somewhat artificially in those regions, 
this time-stepping scheme is no longer effective in simula- 
tions covering a much larger dynamic range. 

In state-of-the-art computer simulations, structures can 
be resolved by A^vir ~ O(IO^) which results in a resolution 
scale of r^es ~ O(10~^ ^vir)- By comparing the dynamical 
scales at rvir and rrcs, we get rdyn(rvir)/rdyn(rrcs) ~ O(IO^) 
or even larger for future high resolution simulations. This 
large dynamical range in time-scales demands an acutely 
adaptive criteria so that the dense, dynamical active regions 
are resolved correctly while the simulations remain fast. 

However, in a general simulation, such as those used for 
cosmological structure formation, it is not straightforward 
to determine the dynamical time of a given particle. This 
depends on the dominant structure responsible for the or- 
bit of a particle which needs to be quickly determined at 
each time-step as the particle is advanced. In this paper we 
develop a new fast method of determining each particle's 
true dynamical time using information directly computed in 
the force evaluation stage of the simulation. This is quite 
different from using the local dynamical time which fails 
dramatically under many circumstances, e.g. consider using 
the local density near the earth to estimate its time-step! 

We present in Section [2] the basic ideas and our im- 
plementation of the new adaptive dynamical time-stepping 
criterion into a tree-code. Detailed and extensive tests are 
presented in Section [S] In Section [4] we conclude. 



structure that the particle is orbiting about. Within colli- 
sionless cosmological simulations this could be some super- 
cluster scale structure, or an individual dark matter halo, 
or some substructure within a dark matter halo. Ideally, we 
would scan the whole sky of the particle and determine the 
structure that gives the dominant contribution to its accel- 
eration. From this dominating structure, we could determine 
the enclosed density and hence find the dynamical time of 
the particle. 

Here we have to distinguish two different regimes. First, 
we have the mean field regime, i.e. particles move in a 
(slowly) varying potential that is determined by the total 
mass distribution. The individual particles are only weakly 
influenced by their direct neighbours and their motions are 
dictated by the sum of more distant particles. This is ensured 
by appropriately softening the short range force thereby 
placing an upper limit on the contribution from an indi- 
vidual particle. In this regime we want the enclosed density 
to be set by the globally dominating structure. The second 
regime is the gravitational scattering regime where we would 
like to follow large angle scattering due to gravitational in- 
teractions, i.e. orbits with eccentricity e ^ 1. Here it is im- 
portant to get the contributions from the closest neighbours 
which dictate the orbital evolution when they are very close 
and when there is little or no force softening. This means 
that the enclosed density is often set by some locally domi- 
nating particle. 

The determination of the enclosed density is quite easy 
for some simple configurations like the 2-body problem or a 
particle orbiting an analytically given spherical symmetric 
structure. However, the generalisation to any given config- 
uration in an A'"-body simulation is not so straightforward 
and we present a simple way in which this can be achieved 
within a tree based gravity code. The specific implementa- 
tion within other code-types may look somewhat different 
but the general scheme and spirit of the method stay the 
same. 



2 DYNAMICAL TIME-STEPPING 

2.1 General idea and description 

In order to advance a particle in a numerical simulation, 
we have to choose a particular time-step for each individual 
particle. Let us consider a particle on a circular orbit in a 
system with spherical symmetric density profile p(r). The 
dynamical time Tdyn('') (or orbital time) of this particle at 
radius r under spherical symmetry is given by 



1 



\/ Gpenc(r) 



where 



Pcnc(r) 



M{r) 



(1) 



(2) 



is the enclosed density within the radius r and M{r) is the 
total mass within radius r. The natural choice for a time- 
step of a particle would therefore be AT cx Tdyn were we not 
faced with the difficulty of determining the enclosed density. 

For a particle in a given landscape of cosmic structure, 
the enclosed density should be set, roughly speaking, by the 



2.2 Implementation within a tree-code 

We use the t ree-code pkdgrav written by Joachim Stadel 
jStadell I2OO1I ') which allows for an adaptive time-stepping 
mechanism where each particle can be on a different time- 
step. The time-steps of the particles are quantized in frac- 
tions of powers of two of a basic time-step To (block time- 
stepping) . Therefore, particles on rung n have an individual 
time-step of 



AT : 



If) 
2" 



(3) 



where Tq is the basic time-step of the simulation and can 
be chosen by the simulator. As stated previously, our time- 
stepping criterion is given by, 

^rD = g<^D / (4) 

where r^D is a free parameter. Therefore, we need to calculate 
the enclosed density ponc for each particle from information 
that is available in a tree-code in order to determine its rung. 

In hierarchical tree-codes, at every time-step two inter- 
action lists are generated for each particle: a list of cells and 
a list of particles that interact with the given particle. The 
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tree structure in such codes is a hierarchical representation 
of the mass in the simulation with each subvolume, or cell, 
being a node in this tree. As we proceed from the root to the 
leaves of this data structure we get an ever finer mass rep- 
resentation of the simulation. The forces from more distant 
cells are calculated by using the multipole expansion of the 
gravitational potential. This expansion makes it clear that 
a finer mass representation, or smaller cell, is required for 
nearby regions than for more distant ones if we want uni- 
form relative errors for each multipole contribution to the 
force. In its simplest form this is realised by a tree-walk al- 
gorithm which, for a given cell, decides whether the use of a 
multipole expansion for this cell satisfies a given error toler- 
ance. If not, this cell is opened and its two or more children 
are considered in the same way. The opening radius of a cell 
which sets an error tolerance is defined by 



(5) 



where 6 is the opening angle and rmax is the distance from 
the centre of mass of the cell to the most distant corner of 
the cell. The numbers are only geometric factors so that in 
the case of a cubic cell with homogeneous density 2/\/3 rmax 
corresponds to the side length of the cube. A cell may only be 
accepted as a multipole interaction if the particle for which 
we are calculating the force is further from the centre of 
mass of the cell than this radius. If a leaf cell (called buckets) 
needs to be opened, then we calculate the interactions with 
each of its particles directly (no multipole expansion is used 
in this case). 

At the end of this procedure each particle in the sim- 
ulation has two interaction lists: 1) a cell list which can be 
thought of as the long range contributions to gravity and 
2) a particle list which accounts for the short range gravi- 
tational interactions. The acceleration and the potential en- 
ergy of each particle are calculated from these two inter- 
action lists. The opening angle varies the ratio of directly 
calculated forces to those calculated via multipole expan- 
sions. It therefore controls force errors and also determines 
the primary cost of a simulation. 

For the calculation of the dynamical time of a particle, 
we generate an additional cell list which provides a reduced 
representation of the particle list, i.e. this list only contains 
the buckets that were opened by the above procedure and 
these buckets are treated in the same way as the cells for the 
long range contributions. The cell list and this additional 
list, which we call the particle-bucket list, form a complete 
tiling of the entire simulation volume except for the local 
bucket of the particle itself, which is not included. 



of which is the dominating region and hence set the cor- 
rect dynamical time-step. The enclosed density for a cell is 
defined by, 

Mc 



Pcnc 



rpc' 



(6) 



where Mc is the total mass of the cell and rpc the vector 
from the location of the particle to the centre of mass of the 
cell. 

(ii) For each of these centres, add up all the pcnc values 
from the other cells in the list that satisfy both of the fol- 
lowing criteria: 



Fpcl 
0.75 



2 IrpCmaxI 
J'PCmax • rpC 



rpcn 



\rpc\ 



(7) 

(8) 



where rpc is the vector from the location of the particle to 
the centre of mass of the cell and rpcmax is the vector from 
the particle to one of the maxima. That means, we add up 
all the Pcnc values of cells that lie within a spherical viewing 
cone of opening angle 2a = 2 arccos(0.75) ~ 83° around a 
maximum cell (Cmax) with the particle (P) being the apex 
extending to 2 |rpcmax I Q See also Fig. \T\ for the geometric 
configuration. If the particle would orbit a perfectly spheri- 
cally symmetric halo at radius r then the dynamical relevant 
mass would lie in the sphere of radius r centred at the ge- 
ometric centre of the halo. Therefore, the angle a is chosen 
so that the volume of the sphere 



Vs 



47r o 



equals the volume of the spherical cone 
Vc = ^{2rf[l~cos{a)] 

resulting in 

Vs_ 1 

Vc 



4[1 



(9) 



(10) 



(11) 



This is reached for cos(a) = 0.75. 

(iii) We now have a summed pcnc of mass contributions 
about each maximum. Only the largest of these, Pcnc.MF, is 
used in determining the dynamical time-step of the particle. 

(iv) Add the local density piocai to the enclosed density 
Pcnc, MP of the particle. We do this in order to account for 
possible contributions from the local bucket of the particle. 

As prefactor we use a valu e of r^p = 0.03 . This choice is 
motivated by earlier studies in IStadell Hool) and is further 
discussed in section 13.51 



2.2.1 Mean field regime algorithm 

We only need to calculate the time-step of a particle when we 
evaluate the force acting on it. The dynamical time of a par- 
ticle is then determined according to the following scheme: 

(i) Pick out the 0.5 percentile highest values of penc from 
both the cell and particle-bucket interaction lists. We regard 
this subset of cells as the centres of dominant contributions 
to the acceleration of the particle, otherwise called maxima. 
Once we have added the contributions of the mass surround- 
ing each of these centres we can make a final determination 



2.2.2 Error terms in leapfrog schemes 

In computer simulations, the system is not evolved by the 
true Hamiltonian but by the approximate Hamiltonian 

Ha^ Ho + AT^Hi + AT^Hi + 0{AT^) . (12) 

^ Adding up the pcnc values shows less scattering in the deter- 
mined time-steps than adding up first the masses of the cells and 
then dividing by the total volume. It also correctly accounts for 
softened contributions to the force from the region close to P since 
the Pcnc contributions are reduced there. 
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Figure 1. Viewing cone for the allowed region of cells to be 
accepted by the time-step criterion. Cmax is the location of the 
maximum cell. We accept all cells that are within the cone of 
opening angle 2a ^ 83° with the particle (P) being the apex 
extending to 2 Irpcmaxl- 



In PKDGRAV, particles are evolved by using a kick-drift-kick 
leapfrog scheme. This means 



Ho 

H2 



Ho 
1 

12 



Hk = H 

{{-ffK, ^^d} , -^d} 



24 



(13) 

{{Hu,Hk},Hk} (14) 



where we have split the true Hamiltonian H into a drift {Hu) 
and kick {Hk) part. A detailed derivation and an expression 
for Hi is given in the Appendix. Therefore, the dominant 
error term is the second order term E2 = AT^H2. 

For a 2-body problem, the Hamiltonian is given by 



2/i 2/ir2 r 



(15) 



Ho 



where jj, 



M1M2 



is the reduced mass and A — GM1M2 



M1+M2 

and where Mi and M2 denote the masses of the two parti- 
cles. The problem is described by the two coordinates r and 
(f and their conjugate momenta 



Pr = fJ,r 

2 ■ r 

Pip — fir (f — L 



(16) 
(17) 



Since the coordinate (p is cyclic, its conjugate momentum is 
an integral of motion, i.e., the angular momentum. 



/iaA(l 



(18) 



is conserved. Here a — z^^, i.e. \a\ is the semimajor axis of 
the ellipse (e < 1) or hyperbola (e > 1) and E is the total 
energy of the orbit. By using a symmetrised time-step. 



Ar = r?D 



G (Ml + M2 




(19) 



we can calculate the higher order error term E2 of the ap- 
proximate Hamiltonian for a 2-body problem and evaluate 
it at pericentrer of the particle's orbit. 



1 (l + 2e) rflA 
24 (1 - e) a 



(20) 



We see that the error depends on eccentricity e of the 
orbit. This allows us to correct for the second order error 
and control the error at pericentre. 



2.2.3 Gravitational scattering regime algorithm 

If the interaction of a particle is in the gravitational scatter- 
ing regime (e.g. it is located close to a super-massive black 
hole), we first determine the mean field value pcnc.MF for 
this particle in exactly the same way as described in section 
12.2.11 However, in order to account for gravitational scatter- 
ing events we need to consider the close particle interactions 
in our determination of the dynamical time in more detail. 
The procedure here is the following: we go through the par- 
ticle interaction list and pick out the highest value of 



/OoncGS = C{e) 



Mp + Ml 
rpp|s 



(21) 



where Mp is the mass of the particle. Mi is the mass of the 
particle in the interaction list, rpp is the particle-particle 
distance and 



C(e) 



1 + 2e 



(22) 



is the additional factor that corrects for eccentricity of the 
orbit. The symmetrisation in pcnc.GS is to cover cases where 
unequal mass particles are involved in the interaction. In 
such cases the heavier particle would be on a much larger 
time-step than the interacting partner, resulting in momen- 
tum conservation problems and unphysical behaviour when 
the mass ratio is large. The eccentricity of two interacting 
particles is given by 



2_EL2 
^2 



(23) 



Hence, for each particle we would like to follow in the gravi- 
tational scattering regime, we have calculated the two values 
of Ponc.MF and Pcnc.GS- The larger of these two is then used. 



3 TIME-STEPPING CRITERION TESTS AND 
BEHAVIOUR 

3.1 General properties 

In order to see how the dynamical time-step criterion works, 
we present the time-step distr ibution in fou r dark matter 
haloe s , so c alled Q:/37-models (iHernauistiil990l : lDehne"nlll993l : 
IZhaolll996l ): with 7 = 1.5,1.0,0.5 & 0.0 where 7 is the 
inner slope of the density profile p(r) cx r~^. The outer 
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Figure 2. We plot the time-step criterion distribution for four profiles with central slope ranging from 7 = 0.0 ... 1.5. The results for 
the dynamical and standard time-stepping criteria are shown - between solid and long dashed lines: dynamical time-step criterion (blue) 
and following the short dashed line: standard time-step criterion (red). Time-steps are in units of the dynamical time at the virial radius 
^dyn(^vir) ~ 12 Gyr. Additionally, we plot the theoretical curves for the standard and dynamical time-stepping criterion we expect for 
a spherical symmetric halo with the given profile. It is evident that the dynamical criterion is much more adaptive with radius than the 
standard criterion. The time-steps taken using the standard criterion remain constant over larger spans in radius than with the dynamical 
time-step criterion. For example, in the NFW profile (7 = 1.0), all particles with a distance smaller than 0.1 r^\^ from the centre are 
on a single constant time-step. Both criteria give equal time-steps at the radius r^q = 4.444 kpc 1.5 X 10~^ r\i^ - independent of the 
density profile. For flat central profiles, the time-step increases below the radius where the acceleration has its maximum! One can also 
see the asymptotic radial behaviour of both schemes in the central region given by equations I I26II respectively I I27I I. We also plot the 
time-step distribution of the dynamical time-step criterion in three bins at ri = r^„, 10~^ rvir,10~'^ rvir of width 0.002 in logarithmic 
scale on the right side of each plot. We can see that the dynamical time-step criterion follows closely a band between rj£, = 0.02. ..0.03 
(long dashed and solid lines). In the two fiat cases (7 = 0.0 and 7 = 0.5), the bin at 10"'^ rvir is close to the resolution limit and the 
distribution is therefore quite broad and noisy. 
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slope is always fi = —3. All haloes have a virial mass of 
Mvir = 10^^ /i-^Mq = 1.429 X 10^^ Mq {h = 0.7) and are 
represented by A'vir ~ 7.5 x 10^ particles within the virial 
radius. This virial mass corresponds to a virial radius of 
ryjr ~ 289 kpc. We fi x the concentration of the NFW profile 
ijNavarro et al.lll996l ) to cnfw = 10 and the concentrations 
of the other profiles were chosen so that the maximum circu- 
lar velocity is reached at the same radius in all haloes. The 
softening of the particles was e = 0.1 kpc ~ 3.5 x 10"** r^ir 
in all haloes. 

We compare our new time-step criterion based on the 
dynamical time with the standard criterion commonly used 
in A'^-body simulations. The standard criterion for selecting 
time-steps in TV-body simulations is based on the accelera- 
tion of the particle. The rung, n, and time-step taken, ATs, 
for the standard criterion is given by, 



(24) 



where e is the softening and a the acceleration of the particle. 
By default, a value 77s = 0.2 is generally used. In spherically 
symmetric systems, the radius Teq where the dynamical and 
standard criteria give the same time-step is given by 



4 



(25) 



independent of the form of the density profile. 

In Fig. [2] we plot the time-step criterion distribution of 
the particles for all haloes as a function of distance from 
the centre of the halo: between solid and long dashed lines 
the values for the dynamical time-step criterion (blue) and 
following the short dashed line the standard time-step cri- 
terion values (red). For a better overview, we only plot 0.1 
per cent of the particles randomly selected from the total 
number of particles in each halo. 

As we can see in Fig. (2] the standard criterion follows 
closely the theoretical curve (short dashed) with \a{r)\ cal- 
culated numerically. The dynamical time-step criterion also 
follows closely a band between rju = 0.02. ..0.03 (long dashed 
and solid lines) for the theoretical curve with pane ~ 
The radius of equal time-steps with the parameters above 



results in roq = 4.444 kpc ^ 1.5 x 10 



also 



On the right side of each plot in Fig. [51 we 
plot the time-step distribution in three bins at ri = 
Tvir, 10~^ Tvir, 10~^ Tvir of width 0.002 iu logarithmic scale. 
For the dynamical time-stepping scheme, most of the parti- 
cles lie in the band between the two curves. Of course, we 
expect a spread in our criterion since our add-up scheme 
does not recognize perfectly the geometry of the surround- 
ing structure. We note however that the add-up scheme is 
generally more conservative for the choice of time-step at a 
given radius than the analytical value. In the two flat cases 
(7 = 0.0 and 7 — 0.5), the bin at 10~^ r^ir is close to the 
resolution limit of the halo and the distribution is there- 
fore quite noisy and relatively broad. Most of the effect due 
to this broader distribution at small radii is absorbed by 
the block time-step scheme's way of discretizing the actual 
time-step taken (see also equation [Jjl ; the time-step value 
provided by the criterion is just an upper limit. 

With a rather conservative choice of 77D = 0.03, we sam- 
ple the orbit of a particle on a circular (tangential) orbit with 
at least 2Ti/rj-D ~ 200 steps. The second curve (short dashed) 



used a value of r/D = 0.02 which corresponds to 2Ti/rj-D ~ 300 
steps per circular orbit. 

The situation is a bit more complicated for particles on 
perfect radial orbits. For a homogeneous sphere, the particle 
will oscillate through the centre of the sphere and describe 
therefore a harmonic oscillator with period T = -y/ 3Tv/{Gp) 
where p is the homogeneous density of the sphere If we 
take this value with pcnc — p we get for a complete radial 
orbit between ^3^/0.03 ^ 100 and ^3^/0.02 « 150 steps 
per oscillation. Since our time-step criterion is very adap- 
tive with radius, the dynamical time will decrease in a steep 
density profile when the particle approaches the centre, so 
that the effective number of steps is even higher. 

The main disadvantage of the standard time-step crite- 
rion (|24|) is the bad adaptivity with radius, i.e. the particles 
are distributed over relatively few rungs. Especially in the 
NFW profile, the particles inside about ten per cent of the 
virial radius are all on the same time-step. For fiatter cen- 
tral profiles with 7 < 1 where p{r) oc r the time-steps 
even increase inside the radius where the acceleration has 
its maximum, in clear contradiction to the behaviour of the 
dynamical time! The asymptotic radial behaviour in the cen- 
tral region of the standard time-stepping criterion is given 
by 



ATs oc 



\a{r)\ 



(7<3) 



(26) 



where 7 is the inner slope of the density profile. In contrast, 
the dynamical time-stepping scheme has the following de- 
pendence 



ATi 



D oc 



Mir] 



oc r 2 (7 < 3) 



(27) 



The standard criterion can only obtain the correct choice of 
time-step in the central regions by shifting the above rela- 
tion, either by choosing a small softening for these particles, 
or by reducing rjs- This automatically leads to an immense 
computational expense due to the overly conservative time- 
steps in the outer parts of the halo or even wrong physical 
behaviour due to the choice of too small softening for the 
physical problem (e.g. undesired scattering of particles). The 
radial scaling of the standard criterion makes it ill suited to 
the study of the centre of galaxies and in other situations 
where a very large dynamic range in density needs to be 
evolved correctly. The dynamical time-stepping technique 
we present is a much more universal approach to choosing 
time-steps in self-gravitating systems, since the basic param- 
eters of the method, such as the angle for adding up mass 
contributions, once determined, are kept fixed for different 
simulations. 



3.2 Elliptic 2-body orbits 

In order to quantify the performance of our adaptive dy- 
namical time-stepping criterion in the scattering regime, we 
performed a series of simulations studying the behaviour of 
high eccentricity 2-body Kepler orbits. After choosing the 



^ Note that iBinnev fc Tremaind lll987l ) defines the dynamical 
time as the time needed by the particle to reach the centre which 
means Tjyn = = ^/3w/{16Gp). 
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-4x10- 



Particle 1 Orbit 1 

Particle 1 Orbit 1000 

Particle 2 Orbit 1 

Particle 2 Orbit 1000 




5x10-' 10° 1.5x10° 2x10° 
X / a 



Figure 3. Top left: run i), top right: run j), bottom left: run k), bottom right: run 1) from Table[T] We plot the orbit from to Tk (orbit 
1) and the orbit from 999 Tk to 1000 Tk (orbit 1000) for each run. All runs used the dynamical time-stepping scheme with eccentricity 
correction. We deliberately did not plot the scaling of the two axes constrained in order to make the small deviations visible. 



masses Mi and M2 of the two bodies, all other quantities 
are fixed, i.e. the orbital time of the Kepler orbit is given by, 



^ / «^(2^)V _ 
^ ^' GKhhh V G(A/i -f M2) 

and the initial total energy is calculated by, 

GM1M2 



Eo = -- 



2a 



(28) 



(29) 



We chose a unit system where Newton's gravitational con- 
stant G = 1 and we fix the orbit in the same way, i.e. the 



semi- major axis is always a = 1. The softening e of the two 
particles was set to 0.1 dpcri in all cases where dpcri = a(l — e) 
is the periapsis distance of the Kepler orbit. Pkdgrav treats 
the forces completely Newtonian if the two particles have a 
distance larger than 2e which is therefore always the case in 
these test runs. Initially, the particles were set in a coordi- 
nate system where the centre of mass is at rest at the origin 
and the two particles were at apoapsis configuration along 
the a;-axis. A short summary of the parameters can be found 
in Table [U 
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Table 1. Summary of elliptic 2-body orbit runs a) - p). 



Run 






^ K 




Rn 

^0 




N 




L^Ulil JJdi C WlLli L Li.ll 


a) 


1 





4.443 




-5 X 10-^ 


2.238 X lO-'J 


256 


D / ec 




b) 


1 


0.9 


4.443 




-5 X 10-1 


-7.278 X 10-1 


2110 


D / ec 


m), n), o) 


c) 


1 


0.99 


4.443 




-5 X 10-1 


6.648 X 10-3 


9394 


D / ec 




d) 


1 


0.999 


4.443 




-5 X 10-1 


2.804 X 10-3 


39143 


D / ec 




c) 


10^ 





1.985 X 10- 


-1 


-5 X 10^ 


2.254 X 10-9 


256 


D / ec 




f) 


10^ 


0.9 


1.985 X 10" 


-1 


-5 X 10^ 


-7.014 X 10-1 


2110 


D / ec 




g) 




0.99 


1.985 X 10- 


-1 


-5 X 10^ 


6.648 X 10-3 


9394 


D / cc 




h) 


10^ 


0.999 


1.985 X 10- 


-1 


-5 X 10^ 


2.804 X 10-3 


39142 


D / ec 




i) 


10^ 





6.283 X 10- 


'3 


-5 X 10^ 


2.242 X 10-9 


256 


D / ec 




j) 


10^ 


0.9 


6.283 X 10 


-3 


—5 X 10^ 


—7.218 X 10 


2110 


D / ec 


P) 


k) 


10'' 


0.99 


6.283 X 10- 


3 


-5 X 10^ 


6.653 X 10-3 


9394 


D / ec 




1) 


10** 


0.999 


6.283 X 10- 


-3 


-5 X 10^ 


2.802 X 10-3 


39142 


D / ec 




m) 


1 


0.9 


4.443 




-5 X 10-1 


1.730 X 10-1 


398 


D / nec 


b) 


n) 


1 


0.9 


4.443 




-5 X 10-1 


-1.263 


316 


S 


b) 


o) 


1 


0.9 


4.443 




-5 X 10-1 


1.895 X 10-3 


2107 


S / rjs = 0.029 


b) 


P) 


lO'' 


0.9 


6.283 X 10- 


3 


-5 X 10^ 


1.913 


3197 


S 


,i) 



The columns are mass ratio M1/M2, eccentricity e, orbital time Tk, initial energy Eg, relative energy change (-E1000 ~ ^o)/l^oli 
number of steps during the first orbit A'^, time-stepping scheme: D - dynamical, S - standard, ec - eccentricity correction, nec - no 
eccentricity correction, run to compare with: b) with m), n) and o) and j) with p). 



We let each run evolve for 1000 Tk {Tk was also the ba- 
sic or longest time-step in the block time-stepping scheme. 
To), measured the total energy -Eiooo after the end of the run 
and calculated the relative energy shift (iSiooo — Eo)/\Eo\. 
These values are also listed in Table [1] From Table [T] we see 
that even for a high eccentricity (e = 0.999) orbit we still 
have relative energy conservation on the level of ~ 10^^ per 
orbit. The general behaviour can also be seen in Fig. [3] where 
we plot the orbit from to Tk and the orbit from 999 Tk to 
1000 Tk for each of the different runs i) - 1). The two orbits 
lie nearly on top of each other except for the e — 0.99 case, 
where the energy gain in the integration was the largest, 
we can see a small deviation. When we compare the runs 
with different mass ratios, we see that the relative energy 
conservation is nearly the same, i.e. the relative energy con- 
servation depends only on the geometry of the orbit. 

In order to illustrate the robustness of our method we 
have performed some further tests. In run m), we switched 
off the eccentricity correction in the symmetrised dynamical 
time-stepping (|2ip . i.e. C(e) = 1. We see that the energy 
gain over 1000 Tk is ~ 240 times larger than in the corre- 
sponding run b) with eccentricity correction. This can also 
be seen visually in Fig.|4] Due to the energy gain, the orbits 
of both particles become wider. 

In run n) we tried to resolve the orbit with the standard 
time-stepping criterion given by equation (|24|l . This crite- 
rion depends on the softening length e of the particle and is 
therefore certainly not ideal in the gravitational scattering 
regime. We've chosen the same value as in run b) where we 
had e = 0.1 dperi = 0.01. From Fig. |4]we see immediately 
when comparing run n) and b) that the standard criterion 
can not capture the dynamics of the orbit. The orbits of the 
two particles become more circular and we get a rotation 
of the whole system. If we had chosen a somewhat larger 
softening, so that it is still smaller than half the periapsis 
distance, it would even look worse since the standard time- 
stepping scheme directly depends on the softening length e 



while the dynamical time-stepping scheme would still per- 
form equally well since it has no such dependence. 

Of course, the dynamical time-stepping criterion with 
eccentricity correction uses a lot more steps per orbit than 
the standard criterion with 77s = 0.2. Therefore we tried in 
run o) a run with an equal number of steps per orbit as run 
b). This is reached for a value i] = 0.029. The energy con- 
servation is still not as good as in the case of the dynamical 
time-stepping with eccentricity correction b) and there is a 
small amount of precession of the periapsis. 

The whole situation becomes even worse when we try to 
resolve 2-body orbits with unequal mass particles. Since the 
standard time-stepping is not symmetric due to the asymme- 
try in acceleration, it is not able to resolve a high mass ratio 
2-body orbit correctly and fails completely. This is shown in 
Fig. [5] Although the light particle makes A'^ = 3197 steps in 
the first orbit (here A'^ in run p) denotes only the number of 
steps of the light particle in Table[T} , the heavy particle takes 
a much larger first step than the light particle. Of course 
when then the light particle approaches, the heavy particle 
gets an immense kick, the total energy becomes positive and 
the whole system drifts apart. 



3.3 Hyperbolic 2-body orbits 

In a similar way, we also tested the new dynamical time- 
stepping scheme for hyperbolic orbits. Initial conditions were 
set up such that the line connecting the two particles en- 
closes an angle of j with the semi-major axis (symmetry 
axis) of the hyperbola. The time for the particle to reach 
the periapsis of the orbit is given by 

Th = (30) 

where r{<j)) describes the angle dependent relative separation 
of the two bodies. The initial conditions used the same unit 
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Run b) / e = 0.9 / M/Mg = 10" 
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Figure 4. Top left: run b), top right: run m), bottom left: run n), bottom right: run o) from Table [T] Run b) shows the dynamical 
time-stepping case with eccentricity correction. The orbit is perfectly followed. In run m), we can nicely see the energy gain visually due 
to the lack of the eccentricity correction in the time-step criterion. Run n) where we tried to resolve a e = 0.9 orbit analogue to run b) 
with the standard time-step criterion. As can be seen, the orbit becomes more circular and the orbital plane starts to rotate which is 
completely unphysical. In run o), we see a run where we used a smaller value of ?7g = 0.029 so that the standard criterion initially makes 
an equal number of steps per orbit as the dynamical time-stepping scheme in run b). This helps a lot but it does still not perform as 
good as the dynamical time-stepping scheme. 



system as the elliptic orbit tests and we again set the soften- 
ing to 0.1 dpcri where dpori = a(e — 1) for hyperbolic orbits. 
A summary of the different runs can be found in Table (2] 

In order to get an integrated effect, we mirrored the 
velocities of the particles after 2 Th and let the runs evolve 
in total for 2000 Th in order to get 1000 pericentre passages. 



In Fig. [6] we plot the most extreme case, run I), with 
mass ratio of 10® and eccentricity e = 1.001. Over 1000 
pericentre passages, there is no visible evolution of the orbit 
and the relative energy change is of order O(10~^) per orbit. 
For the other cases with lower eccentricities and mass ratios, 
the new dynamical time-stepping scheme works optimally 
and we do not show the other orbits here. 
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Table 2. Summary of hyperbolic 2-body orbit runs A) - K). 



Run 


Ml /Mo 




^ H 


^0 


^.^1000 ^0)/ 1-^0 \ 


N 




f OTTiTiPiT'ti Asn + h viin 


A) 


1 


1.1 


2.150 


5 X 10-1 


3.009 X 10-3 


1713 


D / ec 




B) 


1 


1.01 


2.274 X 10-2 


5 X 10-1 


2.119 X 10-3 


4935 


D / ec 




C) 


1 


1.001 


6.709 X 10-4 


5 X 10-1 


-1.298 X 10-2 


15177 


D / ec 




D) 


10^ 


1.1 


9.610 X 10-2 


5 X 10^ 


3.009 X 10-3 


1713 


D / ec 




E) 


10'^ 


1.01 


1.U17 X 10 


5 X 10^ 


2.119 X 10 


4935 


D / ec 




F) 


10^ 


1.001 


2.999 X 10-5 


5 X 10^ 


-1.298 X 10-2 


15178 


D / ec 




G) 


10^ 


1.1 


3.041 X 10-^ 


5 X 10^ 


3.009 X 10-3 


1713 


D / ec 


J), K) 


H) 


lO'^ 


1.01 


3.216 X 10-5 


5 X 10^ 


2.118 X 10-3 


4935 


D / ec 




I) 


lO'^ 


1.001 


9.488 X 10-'' 


5 X 10^ 


-1.298 X 10-2 


15177 


D / ec 




J) 


lO'^ 


1.1 


3.041 X 10-^ 


5 X 10^ 


8.336 X 10-1 


296 


D / nec 


G) 


K) 


10'' 


1.1 


3.041 X 10'^ 


5 X 10^ 


1.360 X 10-1 


333 


S 


G) 



The columns are mass ratio Afi/A/2, eccentricity e, time to reach the pericentre Tjj, initial energy Eq, relative energy change 
(ii'iOOO — Eo)/\Eo\, number of steps during the first orbit N, time-stepping scheme: D - dynamical, S - standard, ec - eccentricity 
correction, nec - no eccentricity correction, run to compare with: G) with J) and K). 



Run p) / e = 0.9 / M/M^ = lO^ 



-5x10 



-1x10 ' 



Particle 1 
Particle 2 Orbit 1 
Particle 2 Orbit 1000 




-1x10 



-5x10- 



-IxlO-s 



-5x10° 
X / a 



-1x10-'^ 10-5 
X / a 

I I I I 



Run I) / e = 1.001 / U^/U^ = 10" 



5x10-3 



CO 



-5x10-3 



-1x10- 



Particle 1 Orbit 1 

Particle 1 Orbit 1000 

Particle 2 Orbit 1 

Particle 2 Orbit 1000 




Figure 5. Run p) where we tried to resolve an e = 0.9, high mass 
ratio orbit with the standard time-stepping criterion. The heavy 
particle gets an immense kick and the total energy becomes posi- 
tive. Thus, the heavy particle drags the light particle with it and 
the whole system drifts apart. The standard criterion therefore 
fails completely to follow the orbit correctly. 

We tried again to resolve high mass ratio orbits with- 
out eccentricity correction (run J) and the standard time- 
stepping scheme (run K) . The results can be seen in Fig. [T] 
If we do not correct for the eccentricity (as in run J), the 
particles gain energy and the orbits become wider. In run K) 
we see the behaviour of the standard time-stepping scheme. 
Due to the large mass ratio, the acceleration of the hght 
particle is quite large and therefore it follows a qualitatively 
correct orbit due to the small time-steps. This is similar to 
the case p) of the elliptic 2-body orbits where the light par- 
ticle describes an elliptical orbit about the massive particle, 
even though this massive particle gets a spurious kick. Once 



Figure 6. The most extreme case, run I), with mass ratio of 10® 
and eccentricity e = 1.001. The dynamical time-stepping scheme 
with eccentricity correction is used. There is no significant evolu- 
tion over 1000 repeated pericentre passages. 

again, the massive particle has a completely incorrect oribit 
wandering around in a very large area due to the spurious 
kicks (compare the scales of the inset plots in runs J and K). 
The lack of momentum conservation in such cases results in 
this contrasting behaviour of the two different mass parti- 
cles. The symmetrization of the time-step criterion restores 
momentum conservation between the two particles involved 
in the gravitational scattering event. 

3.4 Cosmological structure formation 

We also tested the performance of the dynamical time- 
stepping scheme in a cosmological structure formation run. 
For this purpose, we used a fiducial simulation of the Virgo 
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Figure 7. Left: run J), rigiit: run K) from Table|2] In run J) wc see tlie energy gain if the orbit is not followed correctly by not using the 
eccentricity correction in the dynamical time-stepping scheme. On the right, we see run K) where we used the standard time-stepping 
scheme. The heavy particle is wandering around in a much larger area than allowed, as can be seen by comparing the scales of the two 
inset plots. 



cluster (|Ghigna et all Il998l ) in a cosmological framework 
with ~ 1, no cosmological constant i.e. SIa = and 
Ho — 50 km Mpc~^. The simulation cube had a box 
length of L = 100 Mpc and the total mass in the cube 
was Mtot = 6.937 X 10^*^ Mq. The cluster was resolved us- 
ing the standard refinement technique (|Katz fc Whiteiri993l : 
iKatz et al.lll994l ) so that the particle mass in the highest res- 
olution region was 8.604 x 10* Mq and the softening length 
of the lightest particles was e — 5 kpc. The total number 
of particles was 1.314 x 10^. The simulation started at red- 
shift z = 69 and we evolved the cluster to redshift z = 
with three different time-stepping schemes: one dynamical 
and one standard time-stepping run and, for comparison, a 
fixed time-step scheme with 300000 time-steps from z = 69 
to z = (this corresponds to 20000 time-steps down to 
redshift z = 5 in the above described cosmology). With 
this choice, the fixed time-step length corresponds approx- 
imately to the smallest time-step chosen by the dynamical 
criterion during the whole run. Only for a few particles, the 
dynamical scheme did choose smaller steps than this fixed 
time-stepping run. 

The virial radius of the resulting cluster was in all cases 
Tvir ~ 2 Mpc (overdensity ~ 200) and had a final mass of 
A^ciustcr ~ 4.3 X 10" Mq. In Fig. m we plot on the top 
panel the radial density profile for the three runs at redshift 
z = 0. Here pnir) is the radial density of the run with the 
dynamical time-stepping scheme, ps{r) the density profile of 
the run with the standard time-stepping scheme and pF{r) 
the profile of the fixed time-step run. In the lower panel, the 
relative difference (p(r-) — pFij')) / pF^r) is also plotted. The 
softening of the highest resolution particles correspond to 
e = 5 kpc ~ 2.5 x 10~^ rvir- As we can see in Fig. [S] the 




10-3 10-2 10-' 10° 



r / r . 

/ vir 

Figure 8. Density profiles of the three runs of the Virgo cluster 
with the different time-stepping schemes at final redshift 2 = 0: 
Pd ('') is the radial density of the run with the dynamical time- 
stepping scheme, ps(r) the density profile of the run with the 
standard time-stepping scheme and pp(r) the profile of the fixed 
time-step run. The profiles are normalised with respect to the 
critical density Pcrit- On the top panel, the absolute values and 
on the lower panel the relative differences (p(r) — Pf(''))/pf('') 
are plotted. 
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Figure 9. Substructure mass function of the Virgo cluster run at redshift z = 5 (left) and z = (right). There is no substantial difference 
between the runs with different time-stepping schemes visible. 



same radial density profile is obtained for the final cluster 
in this cosmological simulation. 

We also compared the substructure mass function at 
redshifts 2 = and z = 5 for the three runs. For that, we 
used the group finding software skicjf] with a linking length 
of 20 kpc = 4 e and a density and number cut so that only 
structures that are virialised and which are represented by 
at least 100 particles are accepted. In Fig. |9] we plot the 
mass function n{M) (number of substructures of mass M) 
as a function of substructure mass M for output at redshifts 
z = 5 and z = Q. There is no substantial difference between 
the mass functions for the different time-stepping schemes. 

For these low resolution runs we do not expect to see a 
significant difference between the three runs since the scale 
at which the standard scheme begins to take an insufficient 
number of time-steps corresponds approximately to the res- 
olution scale of this simulation. This is just a confirmation 
that the dynamical time-stepping scheme also works for the 
extreme dynamics of a cosmological structure formation run. 

3.5 Dependence on parameters 

3.5.1 Softening length e 

The standard time-step criterion (|24|l depends directly on 
the artificial simulation parameter softening length e. There 
is however no physical basis for this definition. Further- 
more, the functional form of the acceleration in centrally 
flat (7 < 1) haloes is problematic and can lead to nonsen- 
sical time-steps if the resolution is high enough (see Fig. 

Even a simple 2-body problem is not treated properly 
by the standard time-stepping scheme, since the time-steps 

^ http:/ /www-hpcc. astro. washington.edu/tools/skid. html 



depend on acceleration which is not symmetric and again 
there is the meaningless dependence on the softening of the 
particles. 

The dynamical time-stepping scheme only depends in- 
directly on the softening length. If two particles are close 
enough such that their forces are softened, we also use the 
softened values for the penc. In this way the scheme also 
determines an appropriate dynamical time-step when the 
Green's function deviates from the Newtonian 1/r. Further- 
more, the new dynamical time-stepping scheme may be used 
without modification in simulations where the softening is 
set to zero, i.e., where the interactions are never softened. 



3.5.2 Force opening angle 9 

The opening angle 6 determines the weighting of directly 
calculated forces to force contributions coming from the mul- 
tipole expansion. This parameter mainly determines the ac- 
curacy of the force. By including the terms in the particle- 
bucket list, the dynamical scheme does not show a significant 
dependence on the choice of the force opening angle 9. 



3.5.3 Cone viewing angle a 

We normalised the viewing angle a so that the volume of the 
sphere and the cone in Fig. [1] are equal. We also tried larger 
values of cos(q) > 0.75 (i.e. smaller angles) but the resulting 
time-step distribution did not follow the theoretical curves 
as well as for the case of cos(a) — 0.75 (especially close to 
the centre). 
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3.5.4 Number of maxima 

Ideally, one would scan the particle's whole sky for the gravi- 
tationally dominating structure. But this would be computa- 
tionally very expensive. With our choice of the 0.5 percentile 
largest pcnc cells in each of the cell and particle-bucket lists 
we find a good compromise between getting the correct dom- 
inating structure (i.e. low scattering of the enclosed density 
values) and computational speed. Having to consider mul- 
tiple maxima is the main factor which makes the dynam- 
ical time-stepping scheme more expensive than the simple 
schemes used to-date. However, if we loosen the strict geo- 
metrical definition of the viewing cone, then faster schemes 
which rely on the hierarchical tree structure when scanning 
the sky for maxima and their surrounding mass become re- 
alisable. Such algorithmic improvements are being investi- 
gated and will be discussed in future work. 



3.5.5 Prefactor rio 

IStadell l|200ir ) performed stability tests for a leapfrog scheme 
in the drift-kick-drift mode. The result was that 2-body or- 
bits became unstable for choices of r]u ^ 0.1. For these tests, 
the choice of time-steps was also based on the dynamical 
time of the 2-body problem. 

We performed similar tests with the kick-drift-kick 
leapfrog scheme and found that e = 0.9 orbits become un- 
stable in the mean field regime (i.e. without eccentricity cor- 
rection) for choices of rjo too large. Of course by choosing 
a smaller value of j^d, one always gets better precision but 
the computational costs become larger. With the choice of 
rju ~ 0.03 we found a compromise between stability and 
computational costs. 

3.6 Efficiency 

In order to quantify the efficiency of the dynamical time- 
stepping criterion in comparison with the standard criterion, 
we can compare the number of force evaluations for a given 
problem. In a spherically symmetric halo, the number of 
particles in a shell at radius r with thickness dr is given by 



dNp 



A-Kp{r)r^dr 

Mvir/A^vir 



(31) 



The number of time-steps per time interval r for each of 
these dNp particles at radius r is given by 



iVo = — \/Gpcnc(r) = — 

77D 77D 



GM{r) 



(32) 



in the dynamical case. In the case of the standard time- 
stepping scheme, this is instead given by 



iVs 



r 



a{r)\ 



T 



GM{r) 



(33) 



We do not account for the actual block time-stepping scheme 
used in pkdgrav for this numerical estimation. 

The number of force evaluations in that infinitesimal 
thin shell is now simply given by 



dFs = dNp TVs 
respectively 



(34) 




10 » 



Figure 10. Number of force evaluations A'^ per Gyr per radius 
for the dynamical scheme dFpt/dr respectively dF^/dr for the 
standard scheme for an NFW profile dark matter halo. For in- 
finitesimal thin shells at radii larger than roq, the standard time- 
stopping scheme has always a factor (given by equation l|36|l ) 
more force evaluations per Gyr. In an NFW profile, the inner 
slope is 7 = 1 and we have therefore the following asymptotic 
behaviour in the centre: dFpi/dr oc r2 and dF^/dr oc . 



dFu = dNp Np, 



The ratio 
dFs 
dFo 



Re 



Vs 



(35) 



(36) 



shows that above roq defined by equation (|25|l . i.e. the radius 
where both time-stepping criteria give the same value, the 
number of force evaluations at a given radius is always a 
factor Re oc larger. 

In Fig. 1101 we plot the number of force evaluations per 
Gyr per radius dFu /dr and dFs /dr for the NFW profile dark 
matter halo used also in Fig. [2] with iVvir = 7.5 x lO". The 
figure shows the asymptotic behaviour of the curves in the 
central region given by 

—— IX r 2T oc r2 (37) 
dr 

respectively 

^xr^-i''''xr\ (38) 
dr 

In the inner region, i.e at distances smaller than rcq from 
the centre, more force evaluations are done in the dynam- 
ical time-stepping scheme, as expected. Here the standard 
scheme does not give small enough time-steps to follow the 
dynamics. Since in most cases of low resolution simulations 
roq (due to a clever choice of softening or rjs) is around the 
resolution scale, the error is typically small. High resolution 
simulations can however become very slow if one wants to 
resolve the central region correctly with the standard time- 
stepping scheme. 
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Table 3. Number of force evaluations per Gyr for an NFW halo 
with different resolutions A^vir- 





7.5 X 10^ 


7.5 X 10* 


7.5 X lO^'' 




7.612 X 10** 


7.612 X 10^° 


7.612 X 10^2 


Fs 


2.336 X 10^ 


7.388 X 10" 


2.336 X lO^* 




3.069 


9.705 


30.69 



In order to illustrate the efficiency gain, we calculate 
the number of force evaluations per Gyr for three different 
NFW haloes that have the same profile and specifications as 
the one used in Figs [5] and [TU] We only changed the number 
of particles: they are iVvir = 7.5x 10**, 7.5x 10* and 7.5 x 10^° 
particles within the virial radius. For the first halo we chose 
a softening length of ei = 100 pc ~ 3.5 x 10~* Tvir and 
we scaled the softening of the other haloes according to the 
scaling of rimp given by the solution of 

nmp = /i(nmp) , (39) 
where h{r) is the mean particle separation defined by 



h{r) 



Mvir/iVvir 



p(r) 



(40) 



In other words, rimp is the distance of the innermost particle 
to the geometrical centre of the halo and scales as 



1 7=1 



(41) 



resulting in £2 = 10 pc and £3 = 1 pc for the other softenings. 

The number of force evaluations per Gyr (r = 1 Gyr) 
is given by 



Fd = 



dFo respectively Fs = 



dFs 



(42) 



and the numerical results can be found in Table [3] We see 
that Fu scales as Fu oc A^vir whereas -Fs scales approxi- 
mately as Fs oc N^i'^^ in the case of an NFW profile halo. 
This specific scaling of Fs is due to the scaling of the soft- 
ening length e oc rimp resulting in the general scaling of 



Fsociv't^'^-' 



(43) 



For a different scaling of the softening, one gets of course a 
different scaling of _Fs, e.g. e oc N~^y^ results in Fs oc 
independent of 7. This again shows the strong dependence of 
the standard scheme on the softening length e and that the 
the dynamical time-stepping scheme is much more efficient 
than the standard scheme for high resolution simulations. 



4 CONCLUSIONS 

We have developed a physically motivated time-stepping 
scheme that is based on the true dynamical time of the 
particle. We also derive an eccentricity correction for a gen- 
eral leapfrog integration scheme. The combination of these 
schemes allows us to follow quite general dynamical systems 
that may contain a mixture of coUisionless and coUisional 
interacting components. Compared to the standard time- 
stepping scheme used in many A''-body codes it has the fol- 
lowing advantages: 



(i) It does not depend directly on ad-hoc parameters such 
as the softening length e. 

(ii) It gives physically correct time-steps in dark matter 
haloes with arbitrary central cusp slopes. 

(iii) It is faster in high resolution simulations. 

(iv) It allows orbits with eccentricity e ^ 1 to be followed 
correctly. 

(v) It allows us to follow complex dynamical systems 
where scattering events may be important. 

The main conclusion is that one should use a time-step 
criterion that is based on the dynamical time. This scheme 
shows the optimum scaling with the number of particles and 
always gives a physically motivated time-step. 
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APPENDIX A: HAMILTONIAN FORMALISM 

For a physical system that is described by the state z = 
{q,p), where q is the coordinate and p the conjugate mo- 
mentum vector, and this system is evolved under a Hamil- 
tonian H , we can write the formal time evolution (Ha milton 
equations) as l|Saha fc Tremaindll99j : lYoshidalUggi ) 



{z,H} 



dz 

m 

where {, } denote the Poisson brackets defined by 



i 

We can define the operator 



dg dh 
dqi dpi 



dg dh 
dpi dqi 



(Al) 



(A2) 



Hz = {z,H} (A3) 
and write down a formal solution to the time evolution 
z{t) = e'^zo. (A4) 



In computer simulations, the system is evolved by using 
a specific scheme in order to update positions and velocities. 
In PKDGRAV, the following leapfrog scheme is used: during 
a time-step AT, first the velocities are updated (kick mode) 
with a time-step of then the new positions are calculated 
(drift mode) using the new velocities with a time-step of AT 
and finally the velocities are updated to the final values at 
Ar with again a half-step of This is therefore called 
the kick-drift-kick mode and can be described by 



^(Ar) = e~ 



-Hk 



ZO 



(A5) 



where we have split the true Hamiltonian into a drift and a 
kick part 



* http://www.zboxl.org 
^ http://www.zbox2.org 
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H = Ht: 



(A6) 



12 
7 

5760 

7 
1440 
]_ 

~360 
1 

"180 



(A9) 

^ [[Hu,Hk] ,Hk] (AlO) 



24 



By using the Baker-Campbell-Hausdorff series, where we can 
calculate t he higher order terms by an elegant method de- 
veloped by iReinsc E l|200nh . we can write the approximate 
operator Ha under which the system is evolved 

z{AT) = e^^«e^^^°e^^-zo = e^™*zo (A7) 

by 

Ha^Ho + AT^H2 + AT^Hi + 0(/\T^) (A8) 
where 

Ho — Hu + Hk = H 
H2 = — [[Hk,H-d^ ,Hn 

Hi = \ \\\Hu,Hk\,Hk\, Hk\ , Hk] (All) 

Hu,Hk\ ,H-d] , Hk] ,Hk] 

Hd, Hk] , Hk] , [Hu, Hk] ] 

Hk,Hu] ,Hk] ,Hu] , Ho] 

-^[[[HK,Hu]'Hn],[HK,Hn]] 

-^[[[[H^'Ho]^Hu],Hn],Ho] . 

Here, [, ] denote the commutator brackets defined by 

[A,B]^ AB~ BA (A12) 

for non-commutative operators A and B. By using the defi- 
nitions of the operators and the Jacobi identity for Poisson 
brackets, we can calculate the approximate Hamiltonian 

'6\ 



Ha = 
where 

Ho -- 
H2 -- 

Hi 



Ho + AT^H2 + AT^Hi 



0{AT°) (A13) 



Hu+Hk = H (A14) 
^ {{Hk, Hn} > ^^d} - ^ {{^^d, Hk} , Hk} (A15) 



5760 
7 



{{{{Hu,Hk},Hk},Hk},Hk} 
{{{{Hu,Hk},Hu},Hk},Hk} 
{{{Ho,Hk},Hk},{Hb,Hk}} 
{{{{Hk,Hu},Hk},Hu},Hu} 
{{{Hk,Hu},Hu},{Hk,Hu}} 
{{{{Hk,Hu},Hu},Hu},Hu} 



(A16) 



1440 
1 

~360 
1 

"I8O 
1 

~360 
1 

"720 

Note that the replacement of the commutator brackets by 
the Poisson brackets is not trivial. The term H4, is not used 
in our method. However, since H4 is not explicitly given in 
the recent literature we present it here for completeness. 
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